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Abstract. Modern developments in microscopy and image processing 
are revolutionising areas of physics, chemistry, and biology as nanoscale 
objects can be tracked with unprecedented accuracy. However, the price 
paid for having a direct visualisation of a single particle trajectory with 
high temporal and spatial resolution is a consequent lack of statistics. 
This naturally calls for reliable analytical tools which will allow one 
to extract the properties specific to a statistical ensemble from just a 
single trajectory. In this presentation we briefly survey different ana- 
lytical methods currently used to determine the ensemble average dif- 
fusion coefficient from single particle data and then focus specifically 
on weighted least-squares estimators, seeking such weight functions for 
which such estimators will possess ergodic properties. Finally, we ad- 
dress the question of the effects of disorder on such estimators. 



1 Introduction 

Single particle tracking (SPT) can be traced back to the classic studies of Jean Bap- 
tiste Perrin on Brownian motion [T] . With the advent of modern experimental tech- 
niques, recent years witnessed an explosion of different approaches aiming at probing 
physical and biological processes at the level of a single molecule. A SPT experiment 
uses computer-enhanced video microscopy to generate the time series of the position 
ri, T2, rjv at times t\, t^, tjv, of an individual particle trajectory in a 

medium, (see, e.g., [MS])- Properly interpreted, the information drawn from a single, 
or a finite number of trajectories, provides insight into the mechanisms and forces 
that drive or constrain the motion of the particle. Nowadays, single particle tracking 
is extensively used to characterise the microscopic rheological properties of complex 
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media [J, and to probe the active motion of biomolccular motors [5]. In biological 
cells and complex fluids, SPT methods have become instrumental in demonstrating 
deviations from normal Brownian motion of passively moving particles (see, e.g., [5] 

rmspino]). 

The reliability of the information drawn from SPT analysis, obtained at high 
temporal and spatial resolution but at expense of statistical sample size is not always 
clear. Time averaged quantities associated with a given trajectory are usually subject 
to large trajectory-to-trajectory fluctuations even for simple diffusion. For proteins 
transported in the cell cytoplasm or molecules in aqueous environments in general, 
Brownian diffusion is the basic transport mechanism [llj . Notwithstanding that dif- 
fusion is ubiquitous in nature, recent studies have suggested that transport on cell's 
membrane [12], as well as in the cytoplasm [10], may not be limited to pure diffu- 
sion, although the microscopic origin of anomalous diffusion remains unclear. For a 
wide class of anomalous diffusions, described by continuous-time random walks, time- 
averages of certain particle's observables are, by their very nature, themselves random 
variables distinct from their ensemble averages |13U14j . For example, the square dis- 
placement time-averaged along a given trajectory differs from the ensemble averaged 
mean squared displacement [Mil 151116] . By analyzing time- averaged displacements of 
a particular trajectory realization, subdiffusivc motion can actually look normal, al- 
though with strongly differing diffusion coefficients from one trajectory to another [141 
[T5l[T5] and showing substantial ageing effects [T7] . Conflicting results in the identifica- 
tion of the underlying transport mechanisms and their characterization has generated 
in recent years a debate on the most appropriate methodology for the determination 
of the diffusion coefficient from SPT data. 

Even though standard Brownian motion is much better understood than anoma- 
lous diffusion, the analysis of its trajectories is far from being as straightforward 
as one might think, and all the above-mentioned troublesome problems persist. For 
instance, in bounded systems, substantial manifestations of trajectory-to-trajectory 
fluctuations in the first passage time phenomena have been recently revealed [T5UT5] . 

Standard fitting procedures applied to long but finite d-dimensional Brownian 
trajectories unavoidably lead to fluctuating estimates D f of the diffusion coefficient, 
which might be very different from the true ensemble average value D, 

D = 1 * J . 1) 
2dt V ' 

Using different fitting procedures, variations by orders of magnitude have been ob- 
served in SPT measurements of the diffusion coefficient for diffusion of the LacI 
repressor protein along elongated DNA [5U], in the plasma membrane [3] or for dif- 
fusion of a single protein in the cytoplasm and nucleoplasm of mammalian cells |21] . 
The dispersion of Df observed from different single particle trajectories results from 
the rather complex environments in which the measurements are performed. Each 
trajectory will have its own thermal history, particle interactions with different im- 
purities, etc. Moreover, the broad histograms for observed D can also be due to blur 
and localization errors, intrinsic of any experimental measurements, as discussed in 
[251I251I23] . 

The broad dispersion of diffusion constant estimates extracted from SPT analysis 
raises several important questions: Docs an optimal methodology able to determine 
the diffusion coefficient from just one single-particle trajectory exist? If the answer 
to this question is positive, then what is the performance of such methodology to 
the finite length and finite precision of the measured trajectory? Clearly, it is highly 
desirable to have a reliable estimator even for the hypothetical pure cases, such as, 
e.g., unconstrained standard Brownian motion. Such an estimator must possess an 
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crgodic property so that its most probable value should converge to the ensemble 
average and the variance should vanish as the observation time increases. 

This is often not the case and moreover, ergodicity of a given estimator is not 
known a priori and has to be tested for each particular form of the estimator. More- 
over, the knowledge of the distribution of such an estimator could provide a useful 
gauge to disentangle the effects of the medium complexity as opposed to variations in 
the underlying thermal noise driving microscopic diffusion. Recently, much effort has 
been invested in the analysis of this challenging problem and several important re- 
sults have been obtained for the estimators based on the time-averaged mean-square 
displacement 25,26,27], mean maximal excursion |28j . or the maximum likelihood 
approximation [22t|29l[30] . 

In this paper we first review, in Section^ some of the existing statistical method- 
ologies used to extract the diffusion coefficient of single-particle Brownian trajecto- 
ries. In Section [31 we focus on a family of weighted least-squares estimators based on 
single-time averages that we have recently introduced, and in Section |4j we obtain the 
distribution function of the estimate Df and study its ergodic properties. When the 
underlying dynamics is not Brownian, the ergodicity properties of the weighted least- 
squares estimators is not guaranteed. As shown in Section [5j where we present new 
results on the estimation of the diffusion coefficient of a Brownian particle moving on 
a random correlated potential. Section [B] contains our final remarks. 



2 Estimating the diffusion coefficient of single trajectories 

Consider a d-dimensional Brownian process B t with variance Var [B t ] = 2dDt, and 
D as defined in Eq. (fTJ). 

We start this section by considering a simple-minded, rough estimate of Df, defin- 
ing it as the slope of the line connecting the starting and the end-points B t of a given 
trajectory, namely Df = B|/2dt. By definition E [Df] = D. A single trajectory diffu- 
sion coefficient D f so defined is a random variable whose probability density function 
P{Df) is the so-called chi-squared distribution with d degrees of freedom, namely 



1 / d 



where _T(-) is the Gamma-function. Clearly this distribution diverges as Df — > for 
d = 1, P(0) is constant for d = 2, and only for d > 2 the distribution has a bell-shaped 
form with the finite most probable value D*j = (1 — 2/d)D. This means that, e.g., 
for d = 3, the most likely value extracted from a single Brownian trajectory is only 
D f = D/3. 

A more refined method, known as Least-Squares Estimator (LSE), than the previ- 
ous one consists in taking not only the starting and end-points, but the least-squares 
estimate of the trajectory in the full time interval, say t e [0,T], namely 

-Flse = / (B? - l(t)f dt , (3) 
Jo 

where in the simplest case, the dynamical law is taken as l(t) — 2cLDlse^- Minimizing 
Q with respect to -Dlse we obtain the LSE as 

D ^-sdh[ B " tdt - < 4 > 
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Another related method, commonly used in the analysis of SPT, experimental data 
consists in the least-squares fitting of the time-averaged square displacement, also 
called Mean Square Displacement (MSD) [3l ll"l"ll2"T] . For a trajectory followed through 
a time interval T, this is defined as 



1 



-/ 



Fusuit) = — — I (B t+S - B s f ds . (5) 
1 ~ 1 Jo 

At short time lags t — > 0, the time average Fmsd coincides with the ensemble average 
E [Bj] , due to the ergodicity of the diffusion processes. However, due to the finite 
length of experimental trajectories, the MSD analysis is performed over a large 
fraction of intervals t, even when it is clear that the small-i behaviour is often re- 
stricted only to a very small interval compared to the total duration of the trajectory 
|29j . Replacing by the MSD trajectory -Fmsd in Eq. we obtain for a linear 
dynamical law the MSD estimator 

Dmse = J F MS D(t)t dt . (6) 

Calculating the MSD is one of the most popular methods for the analysis of SPT ex- 
perimental data. However, this method presents a number of fundamental limitations, 
leading to non reliable estimations of the diffusion coefficient [24] . Some of these limi- 
tations can be improved by considering weighted MSE, with a time dependent weight 
that takes into account the growth of the variance in time Other limitations 
are rooted to the fact that, differently than the LSE, the MSD estimator is a two- 
time function. This fact renders the MSD estimator particularly fragile with respect 
to the localization errors intrinsic to the experimental acquisition of the data. This 
fragility has been only recently recognized and studied by several authors (see e.g. 
[2"2"Il2"3"Il2~4"] ). With no additional errors, the ensemble average of ©, the average 
of Eq. ([5]) over infinitely different trajectories, coincides with E [Bj] [31]. However, in 
experimental data, the position of the particle at any given time is acquired during 
a finite integration time, yielding a static localization error 8. Since, the localization 
error involves the unknown diffusion coefficient D, the ensemble average of Eq. ([5]) is 
modified as [3"2"] 

E{F MSB (t)]=2(dDt + S 2 ) . (7) 

Additional dynamic errors due to motion blur reduce the previous expression by 
2Dt/3, where r is the time interval at which the position of the particle is recorded 
|33j . Full consideration of these errors and modeling of the motion blur has been 
recently studied in [2"2"ll23j . and in [51] an optimized version of the least-squares 
fitting for the MSD has been proposed. The authors of [53] showed as well that 
the variance of this optimized MSE decays inversely proportional to the length of 
the trajectory, which means that for very long trajectories, the estimated value of the 
diffusion constant becomes trajectory independent. Therefore, this two-time estimator 
is ergo die. 

A conceptually different fitting procedure has been discussed in |29j which amounts 
to maximizing the unconditional probability of observing the whole trajectory B(t), 
assuming that it is drawn from a Brownian process with mean-square displacement 
2dDt. This is the maximum likelihood estimate which takes the value of D that 
maximizes the likelihood of B(t), defined as: 



n(4^)-/ 2 ex P (-^|) (8) 



t=0 
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Minimization of the logarithm of Lt with respect to D yields the Maximum Likelihood 
Estimator (MLE) for the diffusion coefficient of a Brownian trajectory as [29] 

£mle = / dt. (9) 
Jo 1 

This estimator was studied in |29| in one dimension, and it was shown that the MLE 
is superior to those based on the LS unweighted minimization. As a matter of fact, 
the distribution of -Dmle not only appears narrower than the distribution of Z?lse, 
resulting in a smaller dispersion, but also the most probable value of the diffusion 
coefficient appears closer to the ensemble average D [29]. More recently, the same 
conclusions were obtained in |30j for arbitrary dimensions. 



3 Weighted least-squares estimators of the diffusion coefficient 

In a recent paper |34) , we have studied a family of least-squares one-time estimators 
defined as 

u a = ^ f u(t)B 2 t dt, (10) 
where tu(t) is the weighting function of the form 

"(*)=/ (11) 

w {t + i) Q y ' 

a being a tunable exponent, (positive or negative), to - a lag time and A a - the 
normalization constant, appropriately chosen in such a way that E{u Q } = 1, so that 

2dD \J (t +t)«J K } 

Such a normalization permits a direct comparison of the effectiveness of estimators 
corresponding to different values of a. 

The estimator u a in Eq. (|10D minimizes the least-squares functionals with a weight- 
ing function uj(t) = (to + t)~ a . Consider a d-dimensional trajectory B t with t £ [0, T]. 
To estimate the diffusion coefficient Dp from this trajectory, one writes the least- 
squares functional of the form 

F= \[^T (&t-MD F t) 2 dt, (13) 

and seeks to minimize it with respect to the value of Dp, considered as a variational 
parameter. Eq. (|13[) is a generalized weighted least-squares functional, with a cer- 
tain weighting function tu(t), which depending on whether it is a decreasing or an 
increasing function of t, will emphasize the short time or the long time behavior of 
the trajectory B t , respectively. 

After minimization of F with respect to Dp we find that 

^=(?/ T *«W B ?)/(^jf **«(*)) ■ (14) 

Finally, identifying the denominator with the normalization A a and choosing u)(t) = 
(t,Q+t)~ a , we recover our definition ([T0|. Furthermore, note that a = — 1 corresponds 
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to the unweighted LSE and a = 1 corresponds to MLE defined in the previous section. 
Therefore, our generalized estimate (|10[) . possesses a nice property to contain, as 
particular cases, two other commonly used estimators. 

For standard Brownian motion, the lag time to, corresponding to the time at 
which the measurement is started, can be set equal to zero. However, it is useful to 
keep the explicit dependence on to since it is equal to the resolution e at which an 
experimental trajectory is recorded as e = to/T. Moreover, in |30j it was found that 
for anomalous diffusion, or for Brownian motion in presence of disorder to plays a 
significant role. 

In |34j , we have studied the family of estimators (|10[) , and determined a unique 
value of a for which u a is ergodic, so that the single trajectory diffusion coefficient 
Dp — > D &s e = to/T ^ 0. In the next section we briefly sketch the derivation of the 
probability distribution function P{u a ) and discuss its properties. 



4 Distribution of the weighted estimators 

The fundamental characteristic property to derive the probability distribution func- 
tion P(u a ) is the moment-generating function <?(er) of the random variable in Eq. (HOD 
defined as 

$(a) = E [exp (-au a )] , (15) 

where a is positive definite parameter, < a < oo. 

Using the fact that a d-dimcnsional Brownian motion can be decomposed into a 
product of its d one-dimensional components, the generating function can be written 

as 

i-T 



<P(a) = G(*) d 

Following [2"9"II50] . we introduce an auxiliary functional 



cxpf-^ / lo(t)B 2 t (i)cIt 



(16) 



*(x,t)=E* jexp^-^^^r^rU , (17) 

where the expectation is for a Brownian motion starting at x at time t. Clearly, 
G{a) = ^(0,0). This functional satisfies the Feynman-Kac type formula 

V(x, t) = E dB l^(x + dBt,t + dt) (\ - aAa ^ ^dt^ | , (18) 

where dB t is an infinitesimal Brownian increment such that EdB{dBt} = and 
EdB{dBf} = 2Ddt, and E^b denotes averaging with respect to the increment dB t . 
Furthermore, expanding the right-hand-side of the latter equation to second order 
in dB t , linear order in dt and performing averaging, we find that ty(x,t) obeys the 
Schrodinger like equation with a harmonic, time-dependent potential 

^ = -D d -^ + ^x>n*,t), (19) 

subject to boundary condition \P(x, T) = 1 for any x. 

The solution of Eq. (TT9")) , and thus of the generating function (fTB")) . was explicitly 
obtained in [34] for a = 2 and a ^ 2. In the latter case, the moment-generating 
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function, to leading order in e = to/T, is gyben by: 



r(v) — U 2, - 
xij V V xi 



-, -d/2 



-d/2 



for a < 2 



for a > 2 



where i> = 1/(2 — fx), I^(z) is the modified Bessel function [37] and 

d(2 - a) d{a - 2) 

* 1 = ^— and X2= 2(a— I)" 



(20) 
(21) 

(22) 



For the particular case of a = 2, we find that the moment-generating function is 
given for arbitrary e explicitly by 



(8+1) 
2<5e(' 5 - 1 )/ 2 
6-1 



1 



8-1 
6 + 1' 



sh ^2a£<r) 



cosh 



2V2oCo : 



1 + e s ) sinh ( x /2a^cr 



-d/2 



(23) 



First we focus on the variance of the estimator given in (|ll)p. The variance Var(u Q ) 
is obtained by differentiating Eqs. (]20p or (|21|) twice with respect to <r and setting 
a equal to zero. For arbitrary a ^ 2 the variance to leading order in e is then given 
explicitly by 

2 ( (2 - a)/(3 - a), a < 2 , 
Var [«„] = - ^ (24) 
d [ (a-2)/(2a-3), a > 2 . 

The result in the latter equation is depicted in Fig. [T]and shows that, strikingly, the 
variance can be made arbitrarily small in the leading in e order by taking a gradually 
closer to 2, either from above or from below. The slopes at a = 2 + and a = 2~ 
appear to be the same, so that the accuracy of the estimator will be the same when 
approaching a = 2 from above or below. Equation (|24[) . although formally invalid for 
a = 2, also suggests that the estimator in Eq. (Ti"U|) with a = 2 possesses an ergodic 
property. 

A word of caution is now in order. Finite-e corrections to the result in Eq. (|24[) 
are of order of 0(e 2 ~ Q ) for 1 < a < 2, which means that this asymptotic behavior 
can be only attained when e <C cxp (— 1/(2 — a)). In other words, in principle, the 
variance can be made arbitrarily small by choosing a closer to 2, but only at expense 
of increasing either the experimental resolution or the observation time T, which is 
clearly seen in Fig. [TJ 

The solid circles in Fig. Q] correspond to numerical results of random walks on a 
3-dimcnsional lattice and computed P(u a ) using Eq. (|10j) from a large ensemble of 
trajectories, for different values of a and different resolution e. For a < 1.5 or a > 2.5, 
the variance computed numerically is well described by Eq. (|24[) and is independent 
of e (Fig. Near a = 2, corrections due to the finite resolution are noticeable, but 
the numerics clearly show that the variance of the distribution P(u a ) decreases as 
e -> 0. 

We now turn our attention to the distribution function P(u a ), which is obtained 
by inverting the Laplace transform in Eq. (|15p with respect to the parameter a: 



P(u a ) = — / da cxp (au a ) <P(a) . 



(25) 
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Fig. 1. Variance of the distribution P(u a ) for different values of a. The dashed curves 
correspond to Eq. (|24|l . The symbols correspond to the values obtained from numerical 
simulations of 3D random walks for (from light to dark) e = 5 x 10~ 5 ( ), 5 x 10 -6 (•) and 
5 x 1CP 7 (•). The solid squares correspond to the values for Var[u_i] (LSE) and Var [ui] 
(MLE) as indicated by the labels. 



where 7 is a real number chosen in such a way that the contour path of integration 
is in the region of convergence of 'P(cr). Since for a =^ 2 all the poles of the moment- 
generating function lie on the complex plane on the negative real tr-axis |34U35j . we 
can set 7 = in (|25p and find, to leading order in e 



. 1 f°° dz cos (zu a — d(j) a (z)/2) 
P(u a ) = - / 

Jo 



d/4, \ 
Pa [Z) 



where, for a < 2, 

Pa {z) = r 2 (u) (^) 

and the phase 4> is given by 



v-l 



bei^x [ 2 



b a {z) = arctg bcr^-i 2 . / — /bcr^-i 2 . / — 

XiJ V V Xi 



while for a > 2 we have 



and 



4>u{z) = arctg ^ber_, y ^2 
where ber M (x) and bei M (a;) are the Kelvin functions |37j . 



— ) /ber_ 
X2J 



bei 2 



2W- 

X2 



(26) 
(27) 
(28) 

(29) 
(30) 
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a < 2 a > 2 




1 2 3 1 2 3 



U a u a 

Fig. 2. The distribution P(u a ) for different a / 2 in three dimensions. The curves from 
the left to the right (darker to lighter) in the left panel correspond to a = — 1 (LSE), a = 1 
(MLE), ot = 3/2 and a = 1.95 (blue), and in the right panel to a — 5, a = 3, a = 2.5 
and a = 2.05. The symbols in the left panel correspond to numerical simulations of random 
walks in three dimensions. In the inset we show the most probable value it* as a function of 

Q . 

Table 1. Most probable value it* and variance of the normalized estimator u a , for different 
values of a and e — ► 0. 



Estimator 


M* 
"n 


Var(u Q ) 


LSE u_i 


ps 0.44 


0.5 


MLE mi 


ps 0.6 


0.33 


weighted BM Mi. 95 


ps 0.94 


ps 0.032 



In FigsJD we plot P(u a ) in Eq. ([25)1 for three-dimensional systems. Indeed, the 
most probable value m* — > 1 when a — > 2 cither from above or from below. Note 
that for any exponent a 7^ 2, m* is smaller than the average value 1. Nevertheless, 
already for a = 1.95 (or a = 2.05) we get the most probable value u* ps 0.94, which 
outperforms the LSE (m* ps 0.44) and the MLE (u* m 0.6). For a = 1.95 the variance 
Var(it Q ) ps 0.032, which is an order of magnitude less than the variances observed for 
LSE (= 0.5) and the MLE (ps 0.33). Similarly to Fig. [2 finite-resolution corrections 
are negligible for a < 1.5 and a > 2.5, and P{u a ) is well described by Eq. ([20)) . 
For a = 1.95 and finite resolution e = 10~ 7 , we obtain a broader distribution and 
with a smaller u* than the corresponding to Eq. (|20[) for infinite resolution. Note, 
however, that the most probable value of P(ui.9s) that we obtain at finite resolution 
is ps 0.84, which outperforms the LSE and MLE for infinite resolution. As a matter 
of fact, it is evident from Fig. [T]that for any finite resolution, at least < 5 x 10 -5 , 
the variance of the weighted LSE u a outperforms the unweighted LSE and MLE at 
infinite resolution. A comparison between LSE, MLE and our weighted LSE is shown 
in Tabled 
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When a = 2 and e = ta/T small but finite we consider a slightly more general 
form for cj(t): 



w(t) 



2£/t§, for t<t , 



1/t 1 



for t <t<T, 



(31) 



where £ is a tunable amplitude. For such a choice, the moment generating function 
is given explicitly by [34] 



<P((t) 



2(Je (*-i)/2N d / 2 \ ^ 



-d/2 



with 



i ± = (5 ± 1) (ch (V^) ± ^=sh (V^)) , 
where <5 = ^TT4ja and 7 = 2/d(£ + ln(l/e)). Differentiating Eq. ([32]), we find 

4 31n(l/e)-3(l-e) + 2(l-eK + £ 2 



Var[u 2 ] 



3d 



(£ + Ml/e)y 



(32) 
(33) 

(34) 



Now it is required to minimize (|34p with respect to £. We note first that Var(it2) 
is a non-monotonic function of £. However, we find that the optimal value of £ is, for 
arbitrary e 

_ (2 + e)ln(l/e)-3(l-e) 
C " Copt ~ ln(l/ £ ) + 6-1 ■ (35) 

This function and its optimal value are shown in Fig. [3] Therefore, corresponding 
optimized variance is 



Var opt ('U2) = 771 



31n(l/e) - 4 + 5e - e 2 



3d ln(l/e) (ln(l/e) + 1 + 2e) - 3(1 - e) 



(36) 



FromEq. flM]) we find that in 3d Var opt (u 2 ) « 0.144, 0.096, 0.082 for e = 10" 3 ,10" 5 , 10~ 6 , 
respectively. When e — > 0, Var opt (M2) vanishes as 



4 1 

Var op t(u 2 ) ~ 



dln(l/e) ' 



(37) 



Therefore, Var opt (M2) can be made arbitrarily small but at expense of a progressively 
higher resolution. In the limit e — > the distribution converges to a delta-function. 
The estimators with a = 2 are the only, in the family defined by Eq. ([TU]) , that possess 
this crgodic property. 



5 Diffusion in the presence of a random potential 

In this section we consider a Brownian motion in a one-dimensional inhomogeneous 
energy landscape, where the disorder is correlated over a finite length £ c . This model 
gives a simple description of diffusion of a protein along a DNA sequence, for instance, 
where the particle interacts with several neighboring base pairs at a time [36j . The 
total binding energy of the protein is assumed to be a random variable. When the 
particle hops one neighboring base further to the right or to the left, its new energy 
is highly correlated to the value it had before the jump. Slutsky et al. [36] modeled 
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this process as a point-like particle diffusing on a one-dimensional lattice of unit 
spacing with random site energies {Ui}, whose distribution is Gaussian with zero 
mean, variance rj 2 and is correlated in space as ((Ui — Uj) 2 ) = 2r/ 2 [l — exp(— \i— j\/£ c )], 
where £ c is a correlation length. At each time step, the particle located at some site 
i jumps to the left or to the right with probabilities pi oc exp[/3([/j — f/i-i)] and 
qi oc exp[/3(Ui — J7j+i)], respectively, where pi + qt = 1. Diffusion is asymptotically 
normal for any disorder strength rj. Nevertheless, the particle can be trapped in local 
energy minima for long periods of time. During an extended intermediate time regime, 
it is observed that first passage properties fluctuate widely from one sample to another 

m- 

In [30] . we studied the probability distribution function P(u a ) with a = 1 for a 
Brownian particle moving in such a disordered potential. There we found that P(uq) 
is strongly affected by the strength of the disorder and indeed, the presence of disorder 
with short-ranged correlations tends to broaden the distribution of the measured D, 
as it presents an additional source of fluctuations. 

Here we generalise the study in [30] to different values of a. In Fig. [4] we show 
numerical results for the Variance Var [u a ] for a Brownian particle moving on a disor- 
dered potential for different disorder strengths rj and different precision e. As expected, 
for a fixed e, the variance grows in proportion with the disorder strength rj. We do 
not expect that for these dynamics U2 will be an ergodic estimator of the diffusion 
coefficient and as a matter of fact, we cannot say a priori, that there exists a value 
a for which u a is ergodic in this case. Nevertheless, it is interesting to note in Fig. U 
that the variance of u a gets smaller around some value a = a*. Moreover, appar- 
ently a* — > 2 as 77 — > 0. This behaviour is more evident in Fig. [5] where we show the 
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Fig. 4. Variance Var [u a ] for a Brownian particle moving on a short range disordered po- 
tential, as a function of a for different values of precision: e= 10~ 4 ( ), 1(T 5 (•) and 1(T 6 
(•). Each panel corresponds to different disorder strengths as indicated by the labels, and 
/3 = 1. The dashed curve corresponds to the pure Brownian motion (77 = 0) for e = 10 -6 . 



dependence of Var [u a ] on a for a fixed precision e = 10~ 6 and different strengths 
of the disorder rj. As before, Var [u a ] seems to attain a minimum for a value of a 
which moves toward 2 as 77 — > 0. A rough estimation of the optimal value a* from the 
numerical data, shows that the dependence of a* on 77 is consistent with the affinc 
law a* = 2 — r\l 5, as shown in the inset of Fig. [S] 




Fig. 5. Variance Var [u a ] for a Brownian particle moving on a disordered potential, as a 
function of a for e = 1CP 6 and different values of the disorder strength: r\ — 3 ( ), 2 (#), 1 
(•) and 0.5 (•). The dashed curve corresponds to the pure Brownian motion (77 = 0). In the 
inset, a rough estimation of the optimal a* as a function of 77. The dashed curve corresponds 
to a* = 2 — 77/5. 



Will be inserted by the editor 



13 



Finally, to highlight the role of the trajectory-to-trajectory fluctuations, we con- 
sider the probability density function P(co a ) of the random variable 

Ua = (1) 7 (2) ' ( 38 ) 
Ua T Ua 

where and «„' are two identical independent random variables with the same 
distribution P(u a ). The distribution P(w a ), introduced recently in [TS] (see also [2E1 
39,40]), is a robust measure of the effective broadness of P(u a ), which probes the 
likelihood of the event that the diffusion coefficients drawn from two different trajec- 
tories are equal to each other. This characteristic property can be readily obtained 
via an expression |41j 



PM = 7Z 1 V2 j°° uduP{u)p(-^v) . 



(39) 



Hence, P{uj a ) is known once we know P(u a ). 

In [3D], we observed that for r\ w 0.8 the distribution P(io\) undergoes a continu- 
ous shape reversal transition - from a unimodal bell-shaped form to a characteristic 
bimodal M-shape one with the minimum at Wi = 1/2 and two maxima approaching 
and 1 at larger disorder strengths. Therefore, indicating that for r\ > 0.8 sample-to- 
sample fluctuations becomes essential and it is most likely that the diffusion coeffi- 
cients drawn from two different trajectories will be different. Here we have computed 
the distribution P{oj a ) for different values of a, r\ and e. We have found that the 
transition of P{uj a ) described above, is robust with respect to a. Furthermore, we 
have also found that the same transition occurs for fixed strength of the disorder and 
varying a. An example of this transition is shown as a surface plot in Fig. |SJ For 
r] = 3, the transition from bimodal to unimodal P(uj a ) occurs around a ~ 2. This 
means that in the presence of disorder, among the estimators (JTUJ, those for smaller 
a present more trajectory-to-trajectory fluctuations than those with larger a. 



6 Conclusions 

To summarize, we have discussed different statistical estimators of the diffusion coeffi- 
cient of a single (or few) Brownian single-particle trajectories. We have differentiated 
between one-time and two-time estimators. In the case of the former, we have dis- 
cussed a family of weighted least-squares estimators and showed that for standard 
Brownian motion there exist a unique representative of this family ui that possesses 
the ergodic property, meaning that the estimated diffusion coefficient converges to its 
true value in the limit of infinite precision e — > 0. 

Moreover, we have sketched the derivation of the full probability distribution 
function of u a in the leading order in e. For a = 2, we have extracted the optimal value 
of the variance as a function of the precision e, and shown that in the limit e — > 0, 
Var [u e ] vanishes in proportion of 1/ ln(T) . This means that for practical purposes 
the methods based on two-time correlation functions can provide better estimators, 
because the variance of the corresponding estimator decays faster, as 1/T, even in 
the presence of localization errors [22l[23ll24] . 

Furthermore, we studied trajectories of a Brownian particle moving in a random 
potential with short-ranged correlations. Our numerical results suggest that there 
exist an optimal value a = a* for which the variance of u a is minimal. We found that 
a* converges to 2 inversely proportionally to the strength of the disorder 77. Moreover, 
we found that the trajectory-to-trajectory fluctuations increase in proportion to 77 and 
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Fig. 6. Distribution P(u) a ) as a function of the disorder strength a, for fixed disorder 
strength 77 = 3 and e = 10~ 5 . The colour ramp increases from red (small values) to blue 
(large values). 

are inversely proportional to a. This means that a robust estimator of the diffusion 
coefficient for the random motion in disordered potentials is that for which the short 
time dynamics is more efficiently emphasized. The question on the existence of an 
ergodic estimator of the diffusion coefficient for disordered potentials deserves further 
investigation. 

CMM is partially supported by the European Research Council and the Academy of Finland. 



References 

1. Perrin J, C. R. Acad. Sci. 146, (1908) 967; Ann. Chim. Phys. 18, (1909) 5. 

2. Brauchle C, Lamb D C and Michaelis J, Eds., Single particle tracking and single molecule 
energy transfer (Wiley- VCH, Weinheim 2010). 

3. Saxton M J and Jacobson K, Ann. Rev. Biophys. Biomol. Struct. 26, (1977) 373. 

4. Mason T G and Weitz D A, Phys. Rev. Lett. 74, (1995) 1250. 

5. Greenleaf W J, Woodside M T and Block S M, Annu. Rev. Biophys. Biomol. Struct. 36, 
(2007) 171. 

6. Weber S C, Spakowitz A J and Theriot J A, Phys. Rev. Lett. 104, (2010) 238102. 

7. Bronstein I et ai, Phys. Rev. Lett. 103, (2009) 018102. 

8. Seisenberger G et ai, Science 294, (2001) 1929. 

9. Weigel A V, Simon B, Tamkun M M and Krapf D, Proc. Natl. Acad. Sci. USA 108, 
(2011) 6438. 

10. Golding I and Cox E C, Phys. Rev. Lett. 96, (2006) 098102. 

11. Saxton M J, Biophys. J. 72, (1997) 1744. 

12. Pederson T, Nature Cell Biol. 2, (2000) E73-74. 

13. Rebenshtok A and Barkai E, Phys. Rev. Lett. 99, (2007) 210601. 



Will be inserted by the editor 



15 



14. Jeon JHet al., Phys. Rev. Lett. 106, (2011) 048103. 

15. He Y, Burov S, Metzler R and Barkai E, Phys. Rev. Lett. 101, (2008) 058101. 

16. Lubelski A, Sokolov I M and Klafter J, Phys. Rev. Lett. 100, (2008) 250602. 

17. Schulz J, Barkai E and Metzler R, Ageing effects in single particle trajectory averages, 
Phys. Rev. Lett. 110 (2013) 020602. 

18. Meji'a-Monasterio C, Oshanin G and Schehr G, J. Stat. Mech., P06022 (2011). 

19. Mattos T G, Meji'a-Monasterio C, Metzler R and Oshanin G, Phys. Rev. E 86, (2012) 
031143. 

20. Wang Y M, Austin R H and Cox E C, Phys. Rev. Lett. 97, (2006) 048302. 

21. Goulian M and Simon S M, Biophys. J. 79, (2000) 2188. 

22. Berglund A J, Phys. Rev. E 82, (2010) 011917. 

23. X. Michalet, Phys. Rev. E 82, (2010) 041914; 83, (2011) 059904. 

24. X. Michalet and A. J. Berglund, Phys. Rev. E 85, (2012) 061916. 

25. Grebenkov D S, Phys. Rev. E 83, (2011) 061117. 

26. Grebenkov D S, Phys. Rev. E 84, (2011) 031124. 

27. Andreanov A and Grebenkov D S, J. Stat. Mech. (2012) P07001 

28. Tejedor V et al., Biophys J 98, (2010) 1364. 

29. Boyer D and Dean D S, J. Phys. A: Math. Gen. 44, (2011) 335003. 

30. Boyer D, Dean D S, Meji'a-Monasterio C and Oshanin G, Phys. Rev. E 85, (2012) 
031136. 

31. Qian H, Sheetz M P, and Elson E L„ Biophys. J., 60 (1991) 910. 

32. Martin D, Forstner M, and Kas J, Biophys. J., 83 (2002) 2109. 

33. Savin T and Doyle P S, Biophys. J., 88 (2005) 623. 

34. Boyer D, Dean D S, Meji'a-Monasterio C and Oshanin G, Phys. Rev. E 86 (2012) 060101. 

35. Boyer D, Dean D S, Meji'a-Monasterio C and Oshanin G, arXi v: 1301 .43741 

36. M. Slutsky, M. Kardar, and L. A. Mirny, Phys. Rev. E 69, (2004) 061903. 

37. Abramowitz M and Stegun I R, Eds., Handbook of mathematical functions (Dover, New 
York 1972). 

38. Eliazar I, Physica A 356, (2005) 207. 

39. Eliazar I and Sokolov I M, Journal of Physics A: Mathematical and Theoretical 43, 
(2010) 055001. 

40. Eliazar I and Sokolov I M, Physica A 391, (2012) 3043. 

41. Meji'a-Monasterio C, Oshanin G and Schehr G, Phys. Rev. E 84, (2011) 035203. 



